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SUMMARY 

The classical EHL point contact problem is solved using a new 
"system-approach," similar to that introduced by Houpert and Hamrock for the 
line-contact problem. Introducing a body-fitted coordinate system, the troublesome 
free-boundary is transformed to a fixed domain. The Newton-Raphson method can then 
be used to determine the pressure distribution and the cavitation boundary subject to 
the Reynolds boundary condition. This method provides an efficient and rigorous way 
of solving the EHL point contact problem with the aid of a supercomputer and a 
promising method to deal with the transient EHL point contact problem. A typical 
pressure distribution and film thickness profile are presented and the minimum film 
thicknesses are compared with the solution of Hamrock and Dowson. The details of the 
cavitation boundaries for various operating parameters are discussed. 


NOMENCLATURE 

a,b semi-major and semi-minor axes of contact ellipse, m 
E The Young modulus of elasticity, Pa 

u - V 


E' 


equivalent elastic constant, Pa — 


2 

E 7 


(1 




f normal force, N 

G dimensionless material parameter, aE ' 

G dimensionless cavitation boundary 

g cavitation boundary function 

H dimensionless film thickness 
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H dimensionless minimum film thickness 

min 

H dimensionless reference film thickness 

o 

h film thickness, m 

h , minimum film thickness, m 

h Q reference film thickness, m 

k ellipticity parameter 

n normal direction 

P dimensionless pressure 

P Roelands pressure-viscosity constant, l,96e8 Pa 

p pressure. Pa 

p h Hertzian pressure, Pa 

u mean entraining velocity, m/s 

x,y cartesian coordinates 

INTRODUCTION 


In the design of nonconformal contact machine elements, knowledge of elasto- 
hydrodynamic lubrication (EHL) is needed. Since the 1970's, several authors have 
presented their results of the point-contact EHL problem. Among them, the following 
Hamrock and Dowson (H.D.) formula (ref. 1) is widely used in the design of many 
machine elements: 


( **min ) H.C 


3.63U°' 68 G 0,49 W”°* 073 (1 


0 . 60 k . 

e ) 


( 1 ) 


For an EHL solution, a nonlinear integro-dif f erent ial equation must be solved, 
the Reynolds equation and the elasticity equation. The nonlinearities are due to: 

(1) the dependence of the lubricant properties, (viscosity and density), on the 
pressure; (2) the dependence of the film thickness on the pressure; and (3) the free 
boundary at the exit region. Even for the hydrodynamic lubrication, since the free 
boundary is dependent on the pressure distribution, the Reynolds equation has non- 
linear characteristics. It is well known to computational lubrication engineers that 
the numerical treatment of the point-contact EHL problem has inherent difficulties. 
One can see how difficulties arise upon careful consideration of the above mentioned 
nonlinearities. For example, one of the major difficulties is the piezoviscous 
effect. At high loads the viscosity of the fluid can vary by 10 orders of magnitude 
within the conjunction, which caused the pressure spikes and numerical difficulties. 


Another difficulty associated with solving the EHL point-contact case is to 
locate the free boundary where cavitation occurs. In the solution presented by H.D. 
Christopher son' s method (ref. 2) was used together with a Gauss-Seidel iterative 
scheme. The essence of this method is to truncate the negative pressures as they 
occur during iteration and the outlet boundary is located automatically. Oh and 
Rhode (ref. 7) solved the point contact EHL problem using a finite element method and 
Newton's method. But, it has been found that the nonnegativity condition was needed 
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Newton's method. But, it has been found that the nonnegativity condition was needed 
to be checked in each iteration and the discrimination between the continuous film 
region and the cavitated region was troublesome. Though the solution can be obtained 
it is unavoidable that the solution is dependent upon the mesh size distribution near 
the boundary. 

Finally, the large amount of computation time and computer memory space are 
concerns in this calculation. The majority of CPU time is devoted to the calculation 
of the elastic deformation. In general, the Gauss-Seidel iterative method requires 
more than one hundred times of iterations to obtain the converged solution. Further- 
more, to obtain a solution for a given load, one additional loop is required to find 
the reference film thickness. 

This all adds up to the fact that it is very difficult to achieve a stable 
solution at relatively high loads and short CPU times. Recognizing this, Houpert and 
Hamrock (ref. 8) devised an elegant scheme for the line contact case that enabled 
higher load calculations and saved on computational time as well. This scheme was an 
adaptation of Okamura (ref. 9) and became known as the "system approach." Using a 
Newton-Raphson algorithm, the pressures, the integration constant, and the reference 
film thickness are found simultaneously. Here advantage has been taken of the fact 
that the one-dimensional Reynolds equation can be integrated analytically to obtain 
dp/dx and in turn used with the Reynolds' boundary conditions to locate the cavi- 
tation boundary. 

To the author's knowledge, the system-approach has not been successfully applied 
to the point-contact problem. Unlike the line-contact case, the two dimensional 
Reynolds equation can not be integrated analytically. However, a successful 
formulation of the system— approach can nevertheless be accomplished by introducing a 
body-fitted coordinate system and transforming the unknown physical boundary into a 
fixed computational boundary. The unknown boundary function becomes a part of the 
system matrix. In addition, the reference film thickness can be calculated 
simultaneously as was done in the line-contact case. This reduces the number of 
visits to the elastic deformation subroutine substantially. However, as was pointed 
out by Lubrecht et al., (ref. 3), computer memory may be a problem since the Jacobian 
matrix is a full matrix due to the elasticity equation. This problem can be overcome 
by using the block tridiagonal approximation of the system matrix. The matrix 
inversion is accomplished by the Thomas algorithm, and there is no need to store the 
whole Jacobian matrix. Furthermore, the f orce-balance loop can be obviated by 
including it in the system equations and solving simultaneously. 

In this paper, the classical EHL point-contact problem is revisited with a new 
formulation: a free boundary value problem using the system-approach described 

above. The minimum film thicknesses are compared with equation (1) and the details 
of the free boundary are discussed. 


2. ANALYTIC FORMULATION 
2.1 Contact Geometry 

Figure 1 shows the physical model of the elliptical contact, where the x-axis 
represents the rolling or sliding direction, x^ and y fi are the inlet boundaries 
and x * g(y ) is the outlet or cavitation boundary. The ellipticity parameter is 
expressed in terms of the curvature difference (T), the elliptic integral of the 
first (F), and the second kind (S) as 
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2F - S ( 1 + T) 
\| S(1 - T) 


( 2 ) 


where 
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Defining r as R /R , equation (2) can be rewritten as 


y * 


1/2 


<r + 1) - - r 
S 

Therefore, given r, the ellipticity parameter can be calculated iteratively. 

2.2 Governing Equations 

Assuming isothermal conditions and that the lubricant is Newtonian, the steady- 
state Reynolds equation for the point contact problem is 


( 3 ) 


a 

ph 3 dp 

a 

+ 

ph 3 9p 

9x 


ay 



12u 


9(/Ph) 

dx 


(4) 


and, using the parabolic approximation for the geometry, the film thickness is 
expressed by 


x y 2 r r 

h(x,y) - h 0 + + + — -I I 

0 2R 2R 7TE J 


p(x' ,y' )dx' dy ' 


^(x - x' ) 2 + 


(5) 


(y - y' ) 2 


The applied normal force may be balanced by the generated hydrodynamic pressure 
distribution, 

( 6 ) 


f “ I Jn p(x,Y)dx dy 
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Applying the Reynolds boundary condition and symmetry condition at the x-axis, 
the boundary conditions are: 

p = 0 at x = x A ; 0 $ y < y g , 

p = 0 at X A = <X< g(y B ) ; y = y 3 , 


P»0;^f»0 at x-g(y);0<y<y B , 
On 


( 7 ) 


dp . y A 

- 0 at x. <x<g(0);y=0. 

9y 

The viscosity-pressure relation is modeled by the Roelands (11) equation, i.e. 



p a 

*o 

\ 

, P 

Z 

ex P ' 

2 

1 + _ 

PoJ 

-1 


( 8 ) 


a * — (In fi Q +9.67) 
Po 


(9; 


and, the Dowson-Higginson relation (12) is used for the density-pressure relation, 

,9 


0.59x10 + 1 . 34p 

p = p Q P J- n ^ 

0 . 59x10 + p 


( 10 ) 


Letting 


2.3 The Dimensionless Equations 


p hR x x y - g 

P - H - , X = _,Y*_,G*_ 

p h b 2 b a b 


a - ap h , p h 


3 f a - p 

p = JL, p = JL 

27Tab n 0 p 0 


the equations (4) to (6) become 


a 

3x 


pH 3 3p 

i a 
+ — 

/ 

^H 3 3p 


k 2 ay 



x r^ H) 


(ii) 


H(X,Y) -H 0 + i (x 2 + Cl y 2 ) j J n 

* n 


J 


P(X' ,Y' )dX' dY' 

(X - X') 2 + k 2 ( Y - Y') 2 


( 12 ) 


where, X 


12/i n u R 

t 0 m x 


J J fl P(X,Y)dX dY 
k "P h R x 


2n 

T 


fa3 Ph 


c, ■ / c , = 

1 r 2 E'b 


( 13 ) 
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The dimensionless parameters used here can be related to those used by Hamrock and 
Dowson as follows, 


87TUk 


3GWk 

a 


27T 

kH ( 1 + r) 
4Sr 


(14) 


a 

R~ 


6Wk S r 


1/3 


7f 1 + r 

when K = 1 (circular contact), c i and c 2 are 1. 


2.4 Coordinate Transformation 

Introducing the body-fitted coordinate system described in reference 13, 

Y B (X - X A ) 


£ 


the following equations are obtained: 
The Reynolds equation 


G - X. 


r/ = Y 


(15) 


Lj ( P, G , H Q ) 


Y b _9_ 

(G - X A ) 2 9 £ 

£6 ' 9 


\ 

9p 

l 9 

9p 

or 

. « a 

G' 9p 

rr 

NJ 

k 2 Qt i 

l W 

k ja ? 

G - X 3£ 


(16) 


k (G - X, 


.) 8 * 


9p 

l q 


£<g') 2 a 


k 2 (G - X A ) 2 9 £ 


*= a 


^ - x < P H > - 0 

l 9f J g - x A a^ 


where, <3' represents dG/dY. 
The film thickness equation 


H It'D - H 0 ♦ 




+ X, 


c iV 


2c, 


(17) 


+ D(PfG) 

2 


The D represents an integral operator which calculates the elastic deformation of 
two solids in contact resulting from the pressure distribution in the fluid film 
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region (fl). In this paper, the technique presented by Chang (ref. 14) is used. This 
method provides and efficient way of evaluating D without lengthy and complex 
mathematical expressions. Since the coordinate transformation can easily be 
implemented to this method, the details of the algorithm are not presented here. 


The force balance equation is 


L 2 (P,G) - 2 J 0 B J 0 B P <M> 


* d£ dr\ - — = o 

Y„ 3 


(18) 


In the above formulation, L^ is the nonlinear partial differential operator and L 
is integral operator. 


3. NUMERICAL METHOD 
3.1 Spatial Discretization 

To provide a small mesh size near the pressure spike region, an interior 
stretching function (ref. 15) is adopted along the £-axis. The finite difference 
representation of the transformed Reynolds equation is provided in the appendix. 


3.2 Newton's Method 


The system equations are 

L x ( P , G , H Q ) * 0 

L 2 (P,G) - 0 (19 

In reference 13, for hydrodynamic case, has been solved using the Thomas 

algorithm to find P and G, and was used to determine by the force 

balance loop. For the EHL case, the same method can be used. However, if L^ can 
be put in the system equation and be solved simultaneously without creating a 
computer memory problem, the computation will be greatly reduced. 


The system equation for Newton's method can be written as 


where 


f [A] {B}] 

w 



o 

o 

/V 


l L J 

r \ n + 1 

W 

■ « 

+ {5u 


n* 1 

H o 

n 

* H o * 

■ 5h 0 


jA}' 1 - 

2, NI 

- l; 

j - ] 


and [A], {B}, {C} T are the elements in the Jacobian matrix. 


(20) 


( 21 ) 

( 22 ) 


In EHL, due to the integral operator D, [A] is full matrix. But, since large 
amounts of storage and computational time are required to solve it, the block 
tridiagonal matrix approximation can be used. Each block matrix is a full matrix 
which is different from the one-sided arrow shape matrix that resulted from the 
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hydrodynamic case where the elasticity equation is not needed. The unknown G 
matches the residual function at the cavitation boundary where the pressures are 
known from the Reynolds boundary condition. 


Equation (20) can be rearranged as 

[ A ] {<5u} + {b}$h 0 - -(Lj) 

{c) T ($u} . L 2 

From equation (23) 

f<5u} - -[A]" 1 ^} - [A] _1 {b}<5h 0 
Then using equations (24) and (25) 


(23) 

(24) 

(25) 

(26) 

And {<5u> can be calculated using equations (25) and (26). In equation (26), [A] - 1 {B> 
and [A]~ l {L } are obtained by the Thomas algorithm and then stored temporarily and 
used in equation (25). 

The convergence criteria are 

(1) Pressure 




L 2 - {c} [A] 
(c} T [ A] _ 1 (b} 


E E 

i J 


i +1 n 

P - P 

I,J I,J 


< 5.0x10 


-3 


e e 


(2) Cavitation boundary 


Efj 1 - K 


< 5.0x10 


-3 


E° 


a n 
J 


(3) Reference film thickness 


n+ 1 n 

H o - H o 


< 5.0x10 


-3 


H„ 


4. RESULTS AND DISCUSSION 

The dimensionless material parameter used in this analysis is G = 3488 in which 
z = 0.55, p =* 0.018 Ns/m 2 , V = 0.3 and E' = 2.19xl0' u N/m 2 . In figures 2 and 3, 
the pressure distribution and film thickness profile for the circular contact is 
presented. The inlet boundary used for this analysis is defined as X A = -4.0 and 
Y b = 2.0. The maximum pressure is 1.33 times the maximum Hertzian pressure or 
515 MPa which occurs on the x-axis. The dimensionless minimum film thickness is 0.27 
and it occurs at the side-lobes, X = 0.49 and Y = 0.6. 
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The majority of the computation time is used for the calculation of the elastic 
deformations and the differentiations of the residual functions with respect to the 
cavitation boundary function since the integral operator is a function of it. It 
takes about 20 sec on the CRAY-XMP at NASA Lewis Research Center for 1 Newton 
iteration with 3060 nodal points of the whole domain, and, in general, the converged 
solution can be obtained within 3 iterations as long as the initial guess is within 
the sphere of attraction. It was reported (ref. 3) that using the multigrid 
interative method it took 2 hr of CPU time on a VAX11/750 with 2937 nodal points. 

Since a different computer was used, a direct comparison is difficult. The current 
method is quite fast partly because the direct matrix inversion of the block matrices 
is vertorizable which makes it well suited to the supercomputing. Also because the 
amount of visits to the elasticity subroutine is small and there is no need of a 
force balance loop. When the current work is used for transient calculations, the 
previous solution is used as a guess to the next time step and it accelerate the 
solution process, but this is not true for the iterative method. This fact supports 
the current work as a good candidate for transient EHL point contact computation. 

The calculated minimum film thickness in this investigation for various operating 
parameters are provided in table I along with those obtained from the H.D. Formula, 
equation 1. In general, the results from this analysis were higher than those 
predicted by the H.D. for the circular contact case. However, for the elliptical 
contact our results were lower. But the differences do not exceed 10 percent. 

Figure 4 shows another pressure distribution for circular contact where a very 
steep pressure spike occurs. The operating condition is W = 9.154xl0~ 8 , 

U = 1.62xl0“ u , or 0L = 5.723 and \ = 0.862. The maximum pressure is 2.89 times the 
maximum Hertzian pressure or 1.04 GPa. To the authors' experience, the solution is 
so unstable beyond this operation range that the convergence usually fails. When 
U = 6.432xl0‘ 12 , the maximum possible W is 2.367x10’' for circular contact, or 5 = 
7.862 and \ = 0.0964. According to numerous computations, it is found that the 
value of a dictates the numerical stability of current method. The numerical 
stability may be enhanced by reducing the step sizes near the pressure spike region. 
But it should be noticed that the Roelands viscosity-pressure relation is known to be 
valid up to 1 GPa or lower. At such a^high pressure the lubricant behaves as a 
solid-like material and becomes non-Newtonian. Also, recently, it was observed that 
slippage of the lubricant occurs at or very near the surface (ref. 16). Thus it is 
believed that the modification of the classical Reynolds equation including the non- 
Newtonian effect and a more realistic pressure-viscosity relation including the 
thermal effect are needed to investigate the lubrication performance for the high 
load and high speed cases. 

Figure 5 depicts the calculated free boundaries of the circular contacts for 
various operating parameters. The x-axis is somewhat stretched to exaggerate the 
differences. The dotted line represents the Hertzian circle. The general trend is 
that, as expected, the low speed condition results in a curve that conforms more to 
Hertzian (dry) contact circle. Comparing curves 1 and 2, the boundary on the X-axis 
stretches more outside for the higher load but elsewhere it is closer to the Hertzian 
contact circle. Comparing curves 3 and 4, increasing the speed parameter leads to a 
thicker film and tends to straighten out the boundary. 

5. CONCLUSIONS 

The classical EHL point contact problem is solved using a new " system-approach , " 
similar to that introduced by Houpert and Hamrock for the line-contact problem. This 
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requires inverting a system-matrix (i.e., the Jacobian) which via a body-fitted 
coordinate transformation includes boundary conditions at the free boundary. 

Further, a force-balance loop is avoided. Using a Newton-Raphson algorithm, the 
pressures, the cavitation boundary curve, and the reference film thickness are found 
simultaneously. The method is computationally fast and has no problem with locating 
the cavitation boundary. This study revealed that 

1. The minimum film thickness obtained in this study were all within 10 percent 
of the predictions using the H.D. Formula. 

2. The algorithm is well suited to performing transient EHL calculations using 
the supercomputer and the solution at each time step accelerates the succeeding 
solution . 

3. Numerical instabilities were encountered when the value of 5 , that is, W or 
G is high. To obtain a more stable solution, it is believed that the Reynolds 
equation should be modified to include the non-Newtonian effect and a more realistic 
pressure-viscosity relation for high pressure. 

4. The calculated cavitation boundary is near the Hertzian contact circle but 
deviates it for high speed. 
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APPENDIX - DISCRETIZATION ON THE TRANSFORMED REYNOLDS EQUATION FOR k = 1. 
(L 1 ) J = (Rx + £i,i/2 R 7 ) qitl/2 — < P I*1,J " P I.J ) “ ( R 1 + ^I-1/2 R 7 ) < 3 i-1/2,J< P I,J " P I-1,J> 


+ R. 


J+l/2 


( P I(J *1 - P I,j) - " P I'J-1 > 


" R j( R 4 < 3l,J*l/2 DP l ~ R 5^I, 


J-1/2 DP 2 


I- ~ \ 

~ R 6('3l‘l/2,J DP 3 “ ^1-1/2, J DP 4 ) ~ *s\P I* 1 /2 , J H IU/2 , j " Z’l-1/2, J H I-l/2, J ) “ 0 


where, 


2 Y„ 


R i = 


(Gj - X A ) (1 + r^)A£ 

2 


R 2 “ 


(1 + r^A Tj 


R 3 = 


€i 


r^( 1 + r^) (1 + r^AfA rj 
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G J-l/2 “ X A 


Cl G 'j 


(G/ - X A |jr JJ ( I + r € )(l + r v )^T} 


2C I ( G ' J ) ; 


DP, 


(1 + r c )&f (Gj - X A ) 

2Xy b 

R 8 * — 

< G J - X A> A £< 1+r £> 

A C - fi ‘ Ci-x 
« €*.i - Ci 

-Vj - Vj.i 

r A * 7 j*i " V: 

+ R I -i / j+i ) + (rj - 1) (P IfJ + p i,j*i> + < p i*i,o +p i*i,j*i> 
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*I,J 


Pi, J H i,j 
Pi,j 


DP. 


= r g( P i-i,j + P i-l,J-l) + ^ ^ P I,J + P I,J-1^ + ^ P I + 1,J + P I + l,J-l) 


DP, 


" r r/ P I, J-l * P l+l,J-l^ + ^ r Tj ^( P I,J + P I*l,j) + l P I,J^l + P I + 1 # J+1^ 


DP 4 * ” r r? ( p i,j-i * P i-i / j-i^ + ( r r} ^ P I,J + P i-i,J^ + ^ P I,J+1 + P i-l,J*l) 
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TABLE I. - SELECTED MINIMUM FILM THICKNESS 
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e 1 —Point contact EHL model. 
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Figure 2. — Pressure distribution for W = 1.123 x 10 
U = 6.432 x 10~ 12 , G * 3488, and k = 1 . 
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